!--------------------------------------------------------------------------------
! Cumulus Parameterization by G. Grell
!--------------------------------------------------------------------------------

subroutine CUPARM_GRELL_CATT(iens)

  use mem_basic         , only: basic_g
  use mem_tend          , only: tend
  use mem_cuparm        , only: confrq,cuparm_g,cuparm_g_sh

  use node_mod          , only: mynum,   &   ! INTENT(IN)
             		        mxp,     &   ! INTENT(IN)
             		        myp,     &   ! INTENT(IN)
             		        mzp,     &   ! INTENT(IN)
             		        ia,      &   ! INTENT(IN)
             		        iz,      &   ! INTENT(IN)
             		        ja,      &   ! INTENT(IN)
             		        jz,      &   ! INTENT(IN)
             		        i0,      &   ! INTENT(IN)  
             		        j0	     ! INTENT(IN) 
  use mem_grid          , only: time,    &   ! INTENT(IN)
            		   	initial, &   ! INTENT(IN)
            		   	dtlt,	 &   ! INTENT(IN)
            		   	itime1,  &   ! INTENT(IN)
            		   	ngrid,   &   ! INTENT(IN)
            		   	grid_g,  &   ! INTENT(IN)  	  
            		   	dtlongn, &   ! INTENT(IN)
           		   	deltaxn, &   ! INTENT(IN)
           		   	deltayn, &
				npatch
				

  use rconstants        , only: tkmin 
  use extras            , only: extra3d,extra2d 
  use mem_turb          , only: turb_g
  use mem_micro         , only: micro_g
  use mem_scratch       , only: scratch
  use mem_scalar        , only: scalar_g
!srf
  use io_params         , only: frqanl
  use mem_leaf          , only: leaf_g
  use micphys           , only: level

!- use modules for Grell Parameterization
  use mem_grell_param   , only: mgmxp,mgmyp,mgmzp,maxiens,ngrids_cp
  use mem_scratch1_grell, only: ierr4d,jmin4d,kdet4d,k224d,kbcon4d,ktop4d,kpbl4d,kstabi4d, &
                                kstabm4d,xmb4d,edt4d,enup5d,endn5d,deup5d,dedn5d,zup5d,    &
                                zdn5d,iruncon,zcup5d,pcup5d,prup5d,clwup5d,tup5d
  use mem_grell         , only: grell_g,grell_g_sh

!ML -- In case you want to output massflux
  use mem_stilt         , only: imassflx
!ML
  
  implicit none

  integer, intent(IN) :: iens
  integer,parameter :: CPTIME = 0. !orig: CPTIME = 7200.

  integer,parameter :: i_forcing = 1


  !------------------------ deep convection --------------------------------------------
  if(iens == 1) then 
  
  !
  !        Zero out tendencies initially
     if (TIME.eq.0.) then 
        call azero(mxp*myp*mzp,cuparm_g(ngrid)%thsrc(1,1,1))
        call azero(mxp*myp*mzp,cuparm_g(ngrid)%rtsrc(1,1,1))
        call azero(mxp*myp*mzp,extra3d(1,ngrid)%d3(1,1,1)  )
        call azero(mxp*myp    ,extra2d(3,ngrid)%d2(1,1)    )  
     end if 

     !_old : IF(INITIAL.EQ.2.AND.TIME.LT.CPTIME)RETURN
     if(INITIAL.eq.2.and.TIME.lt.CPTIME-dtlt) return
     !_old : IF(MOD(TIME+DTLT+.001,CONFRQ).LE.DTLT.OR.time.lt..01)THEN
     if(mod(TIME,CONFRQ).lt.DTLT.or.time.lt. .01 .or.abs(time-cptime) .lt. 0.01) then !002

        iruncon=1
        
!        print 90,TIME+DTLT,(TIME+DTLT)/3600. +(itime1/100+mod(itime1,100)/60.)
!90      format('   deep conv tendencies updated    time =',F10.1,  &
!             '  real time (hrs) =',F6.1)

        call azero(mxp*myp*mzp,cuparm_g(ngrid)%thsrc(1,1,1))
        call azero(mxp*myp*mzp,cuparm_g(ngrid)%rtsrc(1,1,1))
        call azero(mxp*myp    ,cuparm_g(ngrid)%conprr(1,1) )
        call azero(mxp*myp*mzp,extra3d(1,ngrid)%d3(1,1,1)  )  ! ens output
        call azero(mxp*myp*mzp,extra3d(5,ngrid)%d3(1,1,1)  )  ! cloud/ice tendency
        call azero(mxp*myp    ,extra2d(3,ngrid)%d2(1,1)    )  ! for XFF0 output

        ! >>>> if(mod(TIME-DTLT,frqanl).lt.DTLT) then
        if(mod(time,frqanl) < dtlongn(1) ) then
	     !print*,'---- zerout extra2d array time=',time
             call azero(mxp*myp*mzp,extra3d(4,ngrid)%d3(1,1,1)  )  ! ens output
	endif
	
        
        !srf - use the old way to define the cumulus forcing
	if(i_forcing /= 1) then
             call atob(mxp * myp * mzp,tend%THT(1),grell_g(ngrid)%lsfth(1,1,1))
             call atob(mxp * myp * mzp,tend%RTT(1),grell_g(ngrid)%lsfrt(1,1,1))
        endif	
        call cuparth_catt(                 &
             mynum,                        & !01
             mgmxp,                        & !02
             mgmyp,                        & !03
             mgmzp,                        & !04
             mzp,                          & !05
             mxp,                          & !06
             myp,                          & !07
             ia,                           & !08
             iz,                           & !09
             ja,                           & !10
             jz,                           & !11
             i0,                           & !12
             j0,                           & !13
             maxiens,                      & !15
             iens,                         & !16 
             ngrid,                        & !17
             ngrids_cp,                    & !18
             DTLT,                         & !19
             time,                         & !20
             basic_g(ngrid)%up   (1,1,1),    & !21
             basic_g(ngrid)%vp   (1,1,1),    & !22
             basic_g(ngrid)%wp   (1,1,1),    & !23
             basic_g(ngrid)%theta(1,1,1),    & !24
             basic_g(ngrid)%thp  (1,1,1),    & !24
             basic_g(ngrid)%pp   (1,1,1),    & !25
             basic_g(ngrid)%pi0  (1,1,1),    & !26
             basic_g(ngrid)%dn0  (1,1,1),    & !27
             basic_g(ngrid)%rv   (1,1,1),    & !28
             turb_g(ngrid)%tkep  (1,1,1),    & !29
             tkmin,                          & !30
	     micro_g(ngrid)%rcp(1,1,1),  &! liquid water
    	     micro_g(ngrid)%rrp(1,1,1),  &! pristine
    	     micro_g(ngrid)%rpp(1,1,1),  &
	     micro_g(ngrid)%rsp(1,1,1),  &
    	     micro_g(ngrid)%rap(1,1,1),  &
	     micro_g(ngrid)%rgp(1,1,1),  &
    	     micro_g(ngrid)%rhp(1,1,1),  &
!
	     grid_g(ngrid)%topt  (1,1),      & !29
             grid_g(ngrid)%RTGT  (1,1),      & !30
             !
             grell_g(ngrid)%lsfth(1,1,1)  ,  & !33 
             grell_g(ngrid)%lsfrt(1,1,1)  ,  & !34 
             tend%PT(1),                     & !35
             cuparm_g(ngrid)%THSRC (1,1,1),  & !36 
             cuparm_g(ngrid)%RTSRC (1,1,1),  & !37 
             cuparm_g(ngrid)%CONPRR(1,1),    & !38	
             extra3d(1,ngrid)%d3   (1,1,1),  & !39 ! ensemble output
             extra3d(4,ngrid)%d3   (1,1,1),  & !39 ! ensemble output
             extra3d(5,ngrid)%d3   (1,1,1),  & !39 ! cloud/ice tendency
             extra2d(3,ngrid)%d2   (1,1),    & !39 
             ierr4d,                       & !40
             jmin4d,                       & !41
             kdet4d,                       & !42
             k224d,                        & !43
             kbcon4d,                      & !44
             ktop4d,                       & !45
             kpbl4d,                       & !46
             kstabi4d,                     & !47
             kstabm4d,                     & !48
             xmb4d,                        & !49
             edt4d,                        & !50
             !
             zcup5d,                       & !51 
             pcup5d,                       & !52 
             enup5d,                       & !53
             endn5d,                       & !54
             deup5d,                       & !55 
             dedn5d,                       & !56
             zup5d,                        & !57
             zdn5d ,                       & !58
             prup5d,                       & !59
             clwup5d,                      & !60
             tup5d ,                       & !61
             !
             grell_g(ngrid)%upmf   (1,1),  & !62
             grell_g(ngrid)%dnmf   (1,1),  & !63
             grell_g(ngrid)%xierr  (1,1),  & !64
             grell_g(ngrid)%xktop  (1,1),  & !65
             grell_g(ngrid)%xkbcon (1,1),  & !66
             grell_g(ngrid)%xk22   (1,1),  & !67             
             grell_g(ngrid)%xjmin  (1,1),  & !68
             grell_g(ngrid)%xkdt   (1,1),  & !69
             grell_g(ngrid)%xiact_p(1,1),  & !70
             grell_g(ngrid)%xiact_c(1,1),  & !71
	     confrq,frqanl,                &
	     deltaxn(ngrid)*deltayn(ngrid),&
             leaf_g(ngrid)%patch_area(1,1,1), &
	     npatch,                       &
	     level,                        &
	     grid_g(ngrid)%glat  (1,1),      & !
	     grid_g(ngrid)%glon  (1,1)       ) !

! [ML------------- Stilt - RAMS coupling  ------------------
      if (imassflx == 1) then
        call prep_convflx_to_stilt(mzp,mxp,myp,ia,iz,ja,jz             &
           ,mgmxp,mgmyp,mgmzp,maxiens,ngrid,ngrids_cp                  &    
           ,ierr4d,jmin4d,kdet4d,k224d,kbcon4d,ktop4d,kpbl4d           &
           ,kstabi4d,kstabm4d,xmb4d,edt4d				     &
           ,zcup5d,pcup5d,enup5d,endn5d,deup5d,dedn5d,zup5d,zdn5d      & 
           ,iens)
      end if
! ------------- Stilt - RAMS coupling  ------------------ ML]
     end if 
     call accum(mxp*myp*mzp, tend%tht(1), cuparm_g(ngrid)%thsrc(1,1,1))
     call accum(mxp*myp*mzp, tend%rtt(1), cuparm_g(ngrid)%rtsrc(1,1,1))
     call update(mxp*myp, cuparm_g(ngrid)%aconpr(1,1),cuparm_g(ngrid)%conprr(1,1),dtlt)
     call cupar2mcphysics(mzp,mxp,myp,ia,iz,ja,jz,ngrid,dtlt,& 
                          extra3d(5,ngrid)%d3 (1,1,1),    &
                          basic_g(ngrid)%theta(1,1,1),    & 
                          basic_g(ngrid)%pp   (1,1,1),    & 
                          basic_g(ngrid)%pi0  (1,1,1)     )

  endif 

  !---------------   Shallow cumulus scheme -----------------------------------------
  if(iens == 2) then !006
     if(TIME.eq.0.) then !004
        call azero(mxp*myp*mzp,cuparm_g_sh(ngrid)%thsrc(1,1,1))
        call azero(mxp*myp*mzp,cuparm_g_sh(ngrid)%rtsrc(1,1,1))
     end if 

     if(INITIAL.eq.2.and.TIME.lt.CPTIME-dtlt) return
     if(mod(TIME,CONFRQ).lt.DTLT.or.time.lt. .01 .or. abs(time-cptime).lt. 0.01) then !005
        iruncon=1

!        print 91,TIME+DTLT,(TIME+DTLT)/3600.+(itime1/100+mod(itime1,100)/60.)
!91      format('   shallow conv tendencies updated time =',f10.1,  &
!               ' real time (hrs) =',F6.1)

        call azero(mxp*myp*mzp,cuparm_g_sh(ngrid)%thsrc(1,1,1))
        call azero(mxp*myp*mzp,cuparm_g_sh(ngrid)%rtsrc(1,1,1))           
        call azero(mxp*myp    ,extra2d(2,ngrid)%d2(1,1)       )  


        !srf - use the old way to define the cumulus forcing
	if(i_forcing /= 1) then
             call atob(mxp * myp * mzp,tend%THT(1),grell_g_sh(ngrid)%lsfth(1,1,1))
             call atob(mxp * myp * mzp,tend%RTT(1),grell_g_sh(ngrid)%lsfrt(1,1,1))
        endif	


        call cuparth_shal(       &
             mynum,      &       !01
             mgmxp,      &	 !02
             mgmyp,      &	 !03
             mgmzp,      &	 !04
             mzp,        &	 !05
             mxp,        &	 !06
             myp,        &	 !07
             ia,         &	 !08
             iz,         &	 !09
             !
             ja,         &	 !10
             jz,         &	 !11
             i0,         &	 !12
             j0,         &	 !13
             maxiens,    &	 !15
             iens,       &	 !16
             ngrid,      &	 !17
             ngrids_cp,  &	 !18
             dtlt,       &	 !19
             time,       &       !20
             !
             basic_g(ngrid)%up   (1,1,1),     &   !21
             basic_g(ngrid)%vp   (1,1,1),     &   !22
             basic_g(ngrid)%wp   (1,1,1),     &   !23
             basic_g(ngrid)%theta(1,1,1),     &   !24
             basic_g(ngrid)%thp  (1,1,1),    & !24
             basic_g(ngrid)%pp   (1,1,1),     &   !25
             basic_g(ngrid)%pi0  (1,1,1),     &   !26
             basic_g(ngrid)%dn0  (1,1,1),     &   !27
             basic_g(ngrid)%rv   (1,1,1),     &   !28
             turb_g(ngrid)%tkep  (1,1,1),     &   !29
             !
             tkmin,                          &   !30
	     micro_g(ngrid)%rcp(1,1,1),  &! liquid water
    	     micro_g(ngrid)%rrp(1,1,1),  &! pristine
    	     micro_g(ngrid)%rpp(1,1,1),  &
	     micro_g(ngrid)%rsp(1,1,1),  &
    	     micro_g(ngrid)%rap(1,1,1),  &
	     micro_g(ngrid)%rgp(1,1,1),  &
    	     micro_g(ngrid)%rhp(1,1,1),  &
!
             grid_g(ngrid)%topt     (1,1),   &   !32
             grid_g(ngrid)%RTGT     (1,1),   &   !33
             grell_g_sh(ngrid)%lsfth(1,1,1), &   !34
             grell_g_sh(ngrid)%lsfrt(1,1,1), &   !35
             tend%PT(1),		     &   !36
             cuparm_g_sh(ngrid)%thsrc(1,1,1),&   !37
             cuparm_g_sh(ngrid)%rtsrc(1,1,1),&   !38
             extra3d(2,ngrid)%d3     (1,1,1),&   !39 !<< usando extra3d(2)
             extra2d(2,ngrid)%d2     (1,1),  &   !39 !<< usando extra2d(2)
             !
             ierr4d,                         & !40
             jmin4d,                         & !41
             kdet4d,                         & !42
             k224d,                          & !43
             kbcon4d,                        & !44
             ktop4d,                         & !45
             kpbl4d,                         & !46
             kstabi4d,                       & !47
             kstabm4d,                       & !48
             !
             xmb4d,                          & !49
             edt4d,                          & !50
             zcup5d,                         & !51
             pcup5d,                         & !52
             enup5d,                         & !53
             endn5d,                         & !54
             deup5d,                         & !55
             dedn5d,                         & !56
             zup5d,                          & !57
             zdn5d,                          & !58
             prup5d,                         & !59
             clwup5d,                        & !60
             tup5d,                          & !61
             !
             grell_g_sh(ngrid)%upmf  (1,1),  & !62
             grell_g_sh(ngrid)%xierr (1,1),  & !63
             grell_g_sh(ngrid)%xktop (1,1),  & !64
             grell_g_sh(ngrid)%xkbcon(1,1),  & !65
             grell_g_sh(ngrid)%xk22  (1,1),  & !66
             grell_g   (ngrid)%xierr (1,1),  & !67
	     confrq,frqanl,                  &
	     deltaxn(ngrid)*deltayn(ngrid),  &
             leaf_g(ngrid)%patch_area(1,1,1),&
	     npatch,                         &
	     level)
! [ML------------- Stilt - RAMS coupling  ------------------
       if (imassflx == 1) then
         call prep_convflx_to_stilt(mzp,mxp,myp,ia,iz,ja,jz           &
          ,mgmxp,mgmyp,mgmzp,maxiens,ngrid,ngrids_cp                  &    
          ,ierr4d,jmin4d,kdet4d,k224d,kbcon4d,ktop4d,kpbl4d           &
          ,kstabi4d,kstabm4d,xmb4d,edt4d                              &
          ,zcup5d,pcup5d,enup5d,endn5d,deup5d,dedn5d,zup5d,zdn5d      & 
          ,iens)
       end if
! ------------- Stilt - RAMS coupling  ------------------ ML]

     end if !005
     call accum(mxp*myp*mzp, tend%tht(1), cuparm_g_sh(ngrid)%thsrc(1,1,1))
     call accum(mxp*myp*mzp, tend%rtt(1), cuparm_g_sh(ngrid)%rtsrc(1,1,1))

  end if !006

  !--------- Convective Transport based on mass flux scheme ------------------------------------
  !
  !go to 2000   !desvia do transporte convectivo

  if(iruncon == 1) then

 
     !-------
     !use for testing convective transport (impose no adv and diffuse)
     !if(iens.eq.1) then
     ! call azero(mxp*myp*mzp,A(ISCLT(1))) 
     ! call azero(mxp*myp*mzp,A(ISCLT(2))) 
     ! call azero(mxp*myp*mzp,A(ISCLT(3))) 
     !endif
     !-------
     call azero(mxp*myp*mzp,scratch%scr1(1)) 

     call trans_conv_mflx(iens,scratch%scr1(1))
!-srf -chem
!----- mudar o codigo para salvar tendencias com os arrays da
!----- quimica
     ! Salva tend convectiva 1 no array extra3d
!     if(iens.eq.2) call AE1M1(mxp*myp*mzp,extra3d(3,ngrid)%d3(1,1,1), &
!                              scalar_g(1,ngrid)%sclt(1),              &
!			      extra3d(3,ngrid)%d3(1,1,1)) 
!-srf -chem -end

  end if 
2000 continue

end subroutine CUPARM_GRELL_CATT

!------------------------------------------------------------
subroutine cupar2mcphysics(m1,m2,m3,ia,iz,ja,jz,ngrid,dtlt & 
                          ,sgrell3_3d                &
                          ,theta,pp,pi0)
  use micphys     ,only: level
  use mem_micro   ,only: micro_g
  use rconstants  ,only: cpi
  use mem_scratch
  implicit none
  integer m1,m2,m3,ia,iz,ja,jz,k,i,j,ngrid
  real dtlt
  real, dimension(m1,m2,m3) :: theta, pp, pi0
  real, dimension(m1,m2,m3) :: sgrell3_3d ! liquid/ice tendency from
                                          ! cumulus parameterization

!-srf - nov/2004: not in use
!   return   

   if(level < 2 ) return
   
   if(level == 2) then 
       call mcphysics2(m1,m2,m3,ia,iz,ja,jz,dtlt &
	             ,sgrell3_3d                   &
                     ,micro_g(ngrid)%rcp(1,1,1) )
       return
   endif
   
   if(level == 3)  then 

       call mcphysics3(m1,m2,m3,ia,iz,ja,jz,dtlt,cpi  &
                  , theta, pp, pi0 &
                  ,sgrell3_3d         &        ! cumulus tendency
                  ,scratch%scr1(1)    &        ! tempk
		  ,micro_g(ngrid)%rcp(1,1,1)  &! liquid water
    		  ,micro_g(ngrid)%rrp(1,1,1)  &! pristine
    		  ,micro_g(ngrid)%rpp(1,1,1)  &
		  ,micro_g(ngrid)%rsp(1,1,1)  &
    		  ,micro_g(ngrid)%rap(1,1,1)  &
		  ,micro_g(ngrid)%rgp(1,1,1)  &
    		  ,micro_g(ngrid)%rhp(1,1,1)  &
		  ,micro_g(ngrid)%q6 (1,1,1)  &
    		  ,micro_g(ngrid)%q7 (1,1,1)  )

    endif
return
end  subroutine cupar2mcphysics  

!------------------------------------------------------------------------
subroutine mcphysics2(m1,m2,m3,ia,iz,ja,jz,dtlt,sgrell3_3d,rcp)
implicit none
integer m1,m2,m3,ia,iz,ja,jz,k,i,j
real dtlt
real, dimension(m1,m2,m3) :: rcp,sgrell3_3d

do j = ja,jz
   do i = ia,iz
     do k = 1,m1
       rcp(k,i,j)=rcp(k,i,j)+sgrell3_3d(k,i,j)*dtlt
     enddo
   enddo
enddo
return
end subroutine mcphysics2

!------------------------------------------------------------------------
subroutine mcphysics3(m1,m2,m3,ia,iz,ja,jz,dtlt,cpi,theta,pp,pi0,sgrell3_3d,tempk &
                 ,rcp,rrp,rpp,rsp,rap,rgp,rhp,q6,q7)
implicit none
integer :: m1,m2,m3,ia,iz,ja,jz,k,i,j
real dtlt,cpi
real, dimension(m1,m2,m3) ::  theta,pp,pi0      &
                             ,tempk            &
                             ,rcp,rrp,rpp      &
                             ,rsp,rap,rgp      &
                             ,rhp,q6,q7        

real, dimension(m1,m2,m3) ::  sgrell3_3d


do j = ja,jz
   do i = ia,iz
     do k = 1,m1
            tempk(k,i,j) = theta(k,i,j)*(pp(k,i,j)+pi0(k,i,j))*cpi ! air temp (Kelvin)
    
            if(tempk(k,i,j) .gt. 253.)then
                 rcp(k,i,j)=rcp(k,i,j)+sgrell3_3d(k,i,j)*dtlt
            else
                 rpp(k,i,j)=rpp(k,i,j)+sgrell3_3d(k,i,j)*dtlt
            endif
     enddo
   enddo
enddo

return
end subroutine mcphysics3
!------------------------------------------------------------------------
subroutine prepare_lsf(iwork)

  use mem_grell   ,only: grell_g, grell_g_sh
  use mem_tend    ,only: tend
  use mem_scratch ,only: scratch
  use mem_grid    ,only: time,ngrid,dtlt  
  use mem_cuparm  ,only: confrq ,cuparm_g_sh 
  use node_mod    ,only: mxp,myp,mzp ,ia,iz,ja,jz
  
  implicit none
  character(len=3) :: forcing
  integer,intent(IN) :: iwork

! forcing='RAD' ! you must also coment the lines where lsfcupar is calculated at turb_diff routine
! forcing='PBL' 
! forcing='ADV' 
 forcing='ALL' 


  if(mod(time,confrq).lt.dtlt .or. time .lt. .01) then  

     if(iwork.eq.1) then
        call azero2(mxp * myp * mzp,grell_g_sh(ngrid)%lsfth(1,1,1), grell_g_sh(ngrid)%lsfrt(1,1,1))
        call azero2(mxp * myp * mzp,grell_g   (ngrid)%lsfth(1,1,1), grell_g   (ngrid)%lsfrt(1,1,1))
        return
     endif
     if(iwork.eq.2) then

        call atob(mxp * myp * mzp,tend%THT(1),grell_g_sh(ngrid)%lsfth(1,1,1))! make lsfthsh = tht due radiation
        call atob(mxp * myp * mzp,tend%RTT(1),grell_g_sh(ngrid)%lsfrt(1,1,1))! make lsfrtsh = rtt due radiation
        
	!CALL teste(mzp,mxp,myp,ia,iz,ja,jz,tend%THT(1),grell_g_sh(ngrid)%lsfth(1,1,1),1)
        !CALL teste(mzp,mxp,myp,ia,iz,ja,jz,tend%RTT(1),grell_g_sh(ngrid)%lsfRT(1,1,1),2)

!-srf: for only rad forcing tests
!-srf: you must also coment the lines where lsfcupar is calculated at turb_diff routine
        if    (forcing == 'RAD') then
          call atob(mxp * myp * mzp,grell_g_sh(ngrid)%lsfth(1,1,1),grell_g(ngrid)%lsfth(1,1,1))
          call atob(mxp * myp * mzp,grell_g_sh(ngrid)%lsfrt(1,1,1),grell_g(ngrid)%lsfrt(1,1,1))     
        endif
!- end rad forcing only tests


!-srf: for only turb forcing tests
        if    (forcing == 'PBL') then
         call azero2(mxp * myp * mzp,grell_g_sh(ngrid)%lsfth(1,1,1), grell_g_sh(ngrid)%lsfrt(1,1,1))
         call azero2(mxp * myp * mzp,grell_g   (ngrid)%lsfth(1,1,1), grell_g   (ngrid)%lsfrt(1,1,1))
        endif
!- end turb forcing only tests
!
!
!
        return
     endif
     
!-srf: for rad forcing only tests
!-srf: you must also coment the lines where lsfcupar is calculated at turb_diff routine
     if(forcing == 'RAD' .and. iwork >  2) then
        return
     endif
!- end rad forcing only tests

!-srf: for only turb forcing tests
     if(forcing == 'PBL' .and. iwork ==  3) then
          call atob(mxp * myp * mzp,grell_g_sh(ngrid)%lsfth(1,1,1),grell_g(ngrid)%lsfth(1,1,1))
          call atob(mxp * myp * mzp,grell_g_sh(ngrid)%lsfrt(1,1,1),grell_g(ngrid)%lsfrt(1,1,1)) 
     endif    
     if(forcing == 'PBL' .and. iwork >=  3) then
        return
     endif
!- end turb forcing only tests

     !-srf
     !at this point: PBL term is calculated in diffuse routines and this is the forcing for shallow cumulus
     !CALL teste(mzp,mxp,myp,ia,iz,ja,jz,tend%THT(1),grell_g_sh(ngrid)%lsfth(1,1,1),100)


     ! salva a tendencia total antes da adveccao nos arrais  vt3dm e vt3dn
     ! atencao: certifique -se que os arrays vt3dm e vt3dn nao sao usados(ou zerados) dentro da rotina
     ! de adveccao
     if(iwork.eq.3) then
       
        call azero2(mxp * myp * mzp,scratch%vt3dm(1), scratch%vt3dn(1))
        call   atob(mxp * myp * mzp,tend%THT(1),scratch%vt3dm(1))  ! make vt3dm = tht due all process before adv
        call   atob(mxp * myp * mzp,tend%RTT(1),scratch%vt3dn(1))  ! make vt3dn = rtt due all process before adv

        return
     endif

     if(iwork.eq.4) then
        ! salva a tendencia advectiva nos arrais vt3dm e vt3dn (ae1m1 => a = b - c)
        call ae1m1(mxp * myp * mzp,scratch%vt3dm(1),tend%tht(1),scratch%vt3dm(1)) ! salva tend advectiva thp em a(ivt3dm)
        call ae1m1(mxp * myp * mzp,scratch%vt3dn(1),tend%rtt(1),scratch%vt3dn(1)) ! salva tend advectiva rtp em a(ivt3dn)



!-srf: for only ADV forcing tests
        if(forcing == 'ADV') then
          call atob(mxp * myp * mzp,scratch%vt3dm(1),grell_g(ngrid)%lsfth(1,1,1))
          call atob(mxp * myp * mzp,scratch%vt3dn(1),grell_g(ngrid)%lsfrt(1,1,1))
	  !CALL teste(mzp,mxp,myp,ia,iz,ja,jz,tend%THT(1),scratch%vt3dm(1),3)
          !CALL teste(mzp,mxp,myp,ia,iz,ja,jz,tend%RTT(1),scratch%vt3dn(1),4)
	  return 
        endif    
!- end ADV forcing only tests

	!CALL teste(mzp,mxp,myp,ia,iz,ja,jz,tend%THT(1),scratch%vt3dm(1),3)
        !CALL teste(mzp,mxp,myp,ia,iz,ja,jz,tend%RTT(1),scratch%vt3dn(1),4)


        ! calcula o forcing para conveccao profunda = rad + pbl turb + adv
        ! for deep convection    LSF =  radiation + pbl_turb + advection (ae1p1=> A = B + C)
        call ae1p1(mxp * myp * mzp,grell_g(ngrid)%lsfth(1,1,1),grell_g_sh(ngrid)%lsfth(1,1,1), &
                   scratch%vt3dm(1))
        call ae1p1(mxp * myp * mzp,grell_g(ngrid)%lsfrt(1,1,1),grell_g_sh(ngrid)%lsfrt(1,1,1), &
                   scratch%vt3dn(1))
        
        !CALL teste(mzp,mxp,myp,ia,iz,ja,jz,grell_g(ngrid)%lsfth(1,1,1),grell_g_sh(ngrid)%lsfth(1,1,1),500)
	!CALL teste(mzp,mxp,myp,ia,iz,ja,jz,tend%THT(1),grell_g(ngrid)%lsfrt(1,1,1),6)


        return
     endif


     if(iwork.eq.5) then 

         call accum(mxp*myp*mzp, grell_g(ngrid)%lsfth(1,1,1), cuparm_g_sh(ngrid)%thsrc(1,1,1))
         call accum(mxp*myp*mzp, grell_g(ngrid)%lsfrt(1,1,1), cuparm_g_sh(ngrid)%rtsrc(1,1,1))
     endif
           
  endif

end subroutine prepare_lsf

!---------------------------------------------------

subroutine teste(m1,m2,m3,ia,iz,ja,jz,tht,lsfth,it)

  real, dimension(m1,m2,m3) :: tht,lsfth
  ix=0
!  return

	     print*,'----------------teste:',it
  do j=ja,jz
     do i=ia,iz
	do k=2,m1-1
	   if(j.eq.38.and.i.eq.49)then
	      !IF(tht(k,i,j).gt.1.e-14) THEN
	      !if(k.eq.2) print*,'----------------',it,j,i
	      if(it==1)print*,'TH RAD     =',k,tht(k,i,j)*86400.,lsfth(k,i,j)*86400.
	      if(it==100)print*,'TH RAD+TURB=',k,tht(k,i,j)*86400.,lsfth(k,i,j)*86400.
	      if(it==500)print*,'TH DP SH=',k,tht(k,i,j)*86400.,lsfth(k,i,j)*86400.
	      if(it==2)print*,'RT RAD=',k,tht(k,i,j)*86400.,lsfth(k,i,j)*86400.

	      if(it==30)print*,'TH ADV=',k,tht(k,i,j)*86400.,lsfth(k,i,j)*86400.
	      if(it==40)print*,'RT ADV=',k,tht(k,i,j)*86400.,lsfth(k,i,j)*86400.

	      if(it==3)print*,'TH ADV=',k,lsfth(k,i,j)*86400.
	      if(it==4)print*,'RT ADV=',k,lsfth(k,i,j)*86400.
	      if(it==5)print*,'TH RAD+TURB+ADV=',k,lsfth(k,i,j)*86400.
	      if(it==6)print*,'RT RAD+TURB+ADV=',k,lsfth(k,i,j)*86400.
	   endif
	enddo
     enddo
  enddo
  !IF(ix.eq.10) stop
end subroutine teste

